function [Probm, Ss, rs] = steady_state(Prob, Ss, Hs, rs, tol, max_It);
    global bet0 bet1 tau a b w e pi_theta N_a N_theta N_w N_e w_max w_min alpha...
        xi xi1 ub lb pi_e H w theta Int_pi_w sigmae elasts Pr log_e_mu log_e_sigma fp r_base ss_calc
    fp = 0;
    T = 150;
    ss_calc = 0;
    for t = 1:T
        
        % Find equilibrium spillovers and rental rates
        r_base = 1;
        [Ss, rs] = find_spillovers_ss(Ss, Hs, rs, Prob, tol, max_It);
        
        % Update (w,a) distribution
        Probm = update_dist_ss(Prob, Pr, Ss, rs);        
        my_diff = abs(Probm - Prob);
        eps = sum(my_diff(:));
        
        if eps < 1e-5
            break
        elseif t == T
            error( "Did not converge to SS" )
        end
        Prob = Probm;
  
    end

end